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Abstract:  We  present  an  adaptive  algorithm  which  can  be  used  for  integration 
over  a triangulated  two-dimensional  region  D.  The  integrand  function  may 
depict  some  types  of  singularity  on  subdivision  lines.  The  algorithm  produces 
a sequence  of  approximations  to  the  integral  over  D such  that  an  extrapolation 
to  its  limit  can  be  applied.  The  algorithm  is  a generalization  of  the  TRIEX 
algorithm  by  de  Doncker  and  Robinson. 

Categories  and  subject  descriptors:  G.1.4  Numerical  Analysis:  Quadrature 
and  Numerical  Differentiation  - adaptive  quadrature,  multiple  quadrature;  G.4 
Mathematics  of  Computing:  Mathematical  Software  - algorithm  analysis,  certi- 
fication and  testing;  G.m  Mathematics  of  Computing:  Miscellaneous  - Fortran 
General  Terms,  Theory:  Algorithms 

Additional  Key  Words  and  Phrases:  Automatic  quadrature,  numerical 
multidimensional  integration,  epsilon-algorithm,  singular  integrands 

1.  Introduction 

The  TRIEX  program  [1,2]  was  developed  for  the  automatic  adaptive  quadra- 
ture of  a user-supplied  function  over  a triangle.  The  method  used  in  TRIEX 
is  suited  to  functions  with  some  types  of  singularity  on  the  edges  of  the  given 
triangle  or  on  subdivision  lines  of  the  adaptive  partitioning  strategy.  In  case  of 
a local  difficult  behaviour  of  the  integrand  within  the  triangle,  one  can  proceed 
by  splitting  the  original  triangle  up  so  that  the  problem  occurs  on  boundaries 
of  the  subtriangles  and  call  TRIEX  on  each  of  the  subtriangles  successively. 

However,  setting  the  accuracy  requirements  on  the  separate  triangles  is  hard 
if  a relative  error  has  to  be  satisfied  on  the  entire  region.  Imposing  a too 
restrictive  tolerance  on  triangles  which  have  a small  contribution  to  the  overall 
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integral,  leads  to  unnecessary  work.  Furthermore,  if  the  integral  over  the  entire 
region  is  small  while  the  contributions  over  some  triangles  are  large  in  absolute 
value,  it  may  be  troublesome  to  perform  the  separate  integrations  to  sufficient 
accuracy. 

We  have  extended  and  modified  the  TRIEX  method  so  that  it  will  deal 
with  a triangularized  region,  by  invoking  the  subdivision  process  on  the  set  of 
triangles.  The  resulting  Fortran  subroutine  is  named  TRISET  in  single  and 
DTRIST  in  double  precision.  The  strategy  is  similar  to  that  of  the  quadrature 
routine  (D)QAGP  in  Quadpack  [17],  which  integrates  over  an  interval  that  is 
split  up  at  interior  points  where  some  integrand  difficulties  may  occur.  As  such, 
the  TRISET  algorithm  also  copes  with  various  local  singularities  on  edges  of 
the  given  triangles. 

In  the  present  method,  the  extrapolation  (by  the  e-algorithm  '19])  is  carried 
out  on  a sequence  of  results  each  of  which  is  an  approximation  to  the  integral 
over  the  entire  set  of  triangles  and  corresponds  to  a specific  level  of  subdivision 
of  this  set.  Each  element  in  the  sequence  is  obtained  in  a global  adaptive  way. 
The  algorithm  and  its  theoretical  background  are  outlined  in  the  next  section, 
the  structure  and  use  of  the  program  in  section  3;  some  examples  will  be  given 
in  section  4. 


2.  Algorithm 
Let 

N 

D=  U A*’ 

k-1 

where  the  A*  are  given  non-overlapping  triangles  (by  this  we  mean  that  the 
intersections  of  their  interiors  are  pairwise  disjoint).  Then  an  approximation 

QdU), 

QdU)  ~ JdU)  = / f{x,y)dxdy 
J D 

is  required,  satisfying 

\Io(f)  - QdU)  I < rnax(t0,  tT  |/d(/)|) 

where  ta  and  tT  are  user-prescribed  absolute  and  relative  tolerances  respectively. 
Given  a symmetric  basic  rule  Q for  the  triangle,  we  denote 

q{d]u)  = Y. 

k= i 

where 

for  k = 1, . . . , N 

represents  the  m2-copy  of  Q on  A*,  obtained  by  a scaled  application  of  the  basic 
rule  to  the  subtriangles  in  the  m2-subdivision  of  An  [11,12,13,16].  Thus  each 
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summand  is  a sum  of  integration  results  over  a partition  of  A^  with  subtriangles 
of  size  m~2  times  the  area  of  A*. 

If  each  of  the  quadrature  rule  sums  has  an  error  functional  expansion 

which  allows  an  extrapolation  by  the  e-algorithm  on  the  sequence 

i = 0,1,..., 

then  the  e-algorithm  can  also  be  applied  to  the  sum  expansion  related  to 

{QiP{f)h  i = 0, 1 

Note  that  terms  of  the  same  nature  (dependence  on  vn)  in  the  error  expansions 
arising  from  different  A*,  will  be  eliminated  from  the  sum  expansion  at  the  same 
time.  If,  however,  the  global  expansion  is  complicated,  many  extrapolations  may 
be  needed  to  obtain  the  requested  accuracy. 

A typical  error  expansion  which  allows  the  use  of  the  e-algorithm  for  extrap- 
olation [4]  is  one  with  terms  in 

(log^  Tn)/ma’ , i/j  > 0 integer,  and  ay  > 0 real, 
if  one  extrapolates  on  a sequence  with  a geometric  progression  in  m,  such  as 
m = rrii  — 2l  for  i — 0,1, 

The  asymptotic  expansions  established  for  integration  over  the  simplex  by 
Lyness  [11,12,13]  and  by  Lyness  and  Monegato  [16]  are  of  this  type.  These 
include  two-dimensional  cases  where  /( x,  y)  is  well-behaved  over  the  triangle 
apart  from  possible  singularities  at  the  vertices,  for  example 

f{x,y)  - TpQ(9)h(r)g(x,y) 
over  the  unit  triangle,  with 

r = (z2  + y2)1/2,  9-  tan_1(z/y), 

and  analytic  functions  ©,  h and  g. 

Integrands  with  vertex  and  edge  singularities  of  the  types  studied  in  [14] 
are  also  allowed.  The  latter  papers  deal  with  singularities  of  the  form  xxyprp 
for  integration  over  the  unit  square  (and  the  corresponding  expansions  for  the 
triangle  can  be  derived  from  it).  Edge  singularities  of  the  form  xx  and  xx  logz 
are  handled  in  [18]. 

The  essence  of  the  TRISET  algorithm  is  to  replace  with  a quantity 

that  is  computed  as  a sum  over  fewer  triangles,  not  all  of  the  same  area.  Thus 
TRISET  will  approximate  each  Q ^ \f),  i = 0, 1, . . .,  by  a quantity 

*(/)  = X>L!?(/)  ==  <?£'’(/) 

k = 1 
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where  each  summand  is  a sum  of  integration  results  over  a partition  of  A k with 
subtriangles  of  sizes  > 2-2*  ■ area(Afe)  as  opposed  to  exactly  2-2*  • area(Afc)  in 
the  full  22,-subdivision  of  A*. 

Let  us  refer  to  the  calculation  of  A{(f)  as  “stage  i”  of  the  algorithm.  Stage  0 
simply  involves  the  application  of  Q to  the  given  triangles.  At  the  start  of  each 
further  stage  ( i > 0),  all  triangles  in  the  partition  of  D are  “active”,  meaning 
that  a subtriangle  belonging  to  the  partition  of  Afc  is  of  size  > 2-2t  • area(A*) 
and  hence  can  be  selected  for  subdivision  in  the  course  of  this  stage. 

A stage  (except  for  stage  0)  is  built  up  as  a succession  of  subdivision  steps. 

A subdivision  step  comprises  the  selection  of  the  triangle  with  the  largest  error 
estimate  from  the  active  set,  subdividing  it  into  4 similar  triangles,  computing 
estimates  for  the  integrals  and  errors  over  the  new  subtriangles,  checking  the 
termination  criteria  and  (if  not  terminated)  updating  the  list  of  triangles.  At 
the  end  of  the  i-th  stage,  Q(D  (f)  has  been  successfully  approximated,  i.e.  the 
sum  of  the  error  estimates  over  the  active  set  has  become  small  enough,  and 
A{(f ) is  used  as  a new  extrapolation  entry. 

The  algorithm  can  be  represented  as  given  in  Figure  1.  Details  involving  the 
termination  conditions,  computation  of  the  local  integral  and  error  estimates 
and  the  extrapolation  are  as  in  the  TRIEX  algorithm. 

Algorithm  TRISET ()  { 

Initialize . 

While  global  estimated  error  > global  tolerance  { 

While  estimated  error  over  the  active  set 
> active  set  tolerance  { 

Select  the  next  active  triangle  for  subdivision; 

Subdivide,  and  integrate  over  the  new  subtriangles ; 

Update  the  global  result  and  error  estimate; 

Terminate  in  case  of  abnormal  conditions; 

Update  the  triangle  list; 

} 

Extrapolate  and  estimate  extrapolation  error; 

} 

Terminate  with  either  the  extrapolation  result  or  the  global 
integral  sum  depending  on  their  relative  error  estimates. 

} 


Figure  1:  Algorithm  description 
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3.  Use  and  structure  of  the  program 

The  integrator  is  written  in  Fortran,  and  is  called  as 

CALL  TRISET(F , N , X ,Y .EPSABS .EPSREL .MEVALS , ICLOSE , MAXTRI , 

* RESULT , ABSERR , MEVALS , NTR , IER.WRK , IWRK ) 

TRISET  returns  an  approximation  RESULT  to  the  sum  / of  the  integrals  of  the 
function  F(X,Y)  over  the  given  set  of  triangles,  and  tries  to  satisfy  an  accuracy 
requirement  of  the  form 

\I  - RESULT]  < max{EPSABS,  EPSREL  ■ \I\} 

where  EPSABS  and  EPSREL  are  the  requested  absolute  and  relative  accuracies. 
The  user  specifies  N triangles  by  supplying  the  abscissas  and  ordinates  of  their 
vertices  in  X(3,*)  and  Y(3,*),  respectively.  MEVALS  is  a user-set  limit  on 
the  number  of  F(X,Y)  evaluations.  The  integrator  returns  an  estimate  (often 
pessimistic)  of  the  absolute  error  1 1 — RESU LT\  in  ABSERR.  If  all  goes  well, 

\I  - RESULT]  < ABSERR  < ma x{EPSABS,  EPSREL  ■ |JJ}. 


Figure  2:  Structure  of  TRISET 
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NEVALS  is  the  number  of  function  evaluations  used.  When  problem  situations 
are  detected,  an  error  code  is  issued  in  IER. 

Two  local  quadrature  modules  are  provided  with  the  subroutine.  The  user 
selects  module  LQO  by  setting  the  input  parameter  ICLOSE  to  0.  LQl  is  selected 
by  setting  ICLOSE  to  1 or  2.  Each  local  quadrature  module  is  based  on  Lyness 
and  Jespersen  quadrature  rules  for  the  triangle  [15].  LQl  uses  the  Lyness  and 
Jespersen  rule  of  degree  9 with  19  points  and  the  rule  of  degree  11  with  28 
points.  LQO  uses  the  Lyness  and  Jespersen  rule  of  degree  6 with  12  points 
and  the  rule  of  degree  8 with  16  points.  LQl  is  usually  more  accurate  than 
LQO  for  fairly  well-behaved  functions.  However,  LQO,  unlike  LQl,  uses  function 
values  only  at  interior  points  of  the  triangle,  so  LQO  can  be  used  in  cases  where 
the  integrand  has  singularities  on  an  edge  without  explicitly  defining  F(X,Y) 
there.  When  a closed  rule  is  used  in  this  case,  it  is  common  practice  to  set  the 
function  to  0 where  the  problem  occurs.  In  some  cases  of  vertex  singularities, 
the  related  asymptotic  error  expansions  are  known  to  be  of  the  same  form  as 
the  ones  obtained  with  an  open  rule  [10]. 

With  ICLOSE  = 1,  the  extrapolation  is  performed  as  in  TRIEX,  on  the 
basis  of  the  error  expansion  for  the  integral  contributions  of  degree  9.  This 
alleviates  problems  with  possible  corruption  of  the  extrapolation  entries  in  case 
of  singularities  on  edges,  where  degree-11  rule  evaluations  are  needed.  Note 
that  the  degree-11  rule  does  not  require  evaluations  at  the  vertices.  It  can  be 
observed  that  use  of  the  degree- 11  quadrature  sums  for  the  extrapolation  is 
often  more  efficient  than  the  degree-9  rule  in  the  presence  of,  for  example,  a 
vertex  singularity  or  some  types  of  oscillatory  behaviour  of  the  integrand,  as  is 
explained  by  the  higher  degree  of  the  rule.  ICLOSE  = 2 provides  the  option 
of  extrapolating  on  the  sequence  of  global  quadrature  sums  obtained  with  the 
degree-11  rule.  At  edge  singularities  (of  the  original  triangles  or  subtriangles) 
it  is  assumed  that  the  integrand  is  set  to  zero  by  the  user. 

The  open  rule  used  when  ICLOSE  = 0,  with  less  points  per  rule  evaluation, 
is  in  general  efficient  when  many  subdivisions  are  required  to  construct  each 
extrapolation  entry. 

Refer  to  Figure  2 for  a diagram  of  the  routines  called  by  TRISET.  Module 
TRISET  allots  array  storage  space  and  calls  TRIS1  which  performs  the  work. 
RlMACH  is  the  Bell  Labs  [3]  routine  which  supplies  the  machine  constants 
needed.  XERROR  [7]  is  used  for  error  handling.  TRISET  is  structured  sim- 
ilarly to  TRIEX,  with  the  main  differences  being  the  addition  in  TRISET  of 
TRIS1  and  of  two  local  quadrature  modules  (TRIEX  uses  only  one  closed  rule 
corresponding  to  LQl). 

Further  information  involving  the  usage  of  the  program  is  given  in  the  Pro- 
logue of  the  routine. 

4.  Examples 

In  Tables  2 and  3 we  show  the  results  from  DTRIST  for  the  integration  of 
the  functions  exp(  — x — y)((l  — z )2  + y2)  09  and  y|y(0.5  — z)|(0.5  - z)  respec- 
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tively.  over  the  star-shaped  region  composed  by  the  triangles  listed  in  Table  1 
and  shown  in  Figure  3. 


Triangle 

Vertices 

1 

(-1/3,0)  (1,0)  (0,.5) 

2 

(0,-5)  (1,0)  (.7, .8) 

3 

(.7, .8)  (1,0)  (3tt,1) 

4 

(0,-5)  (.7, .8)  (.6,22) 

5 

(— e,1.7)  (-1/3,0)  (0,.5) 

6 

(-10,-10)  (1/3, -.5)  (-1/3,0) 

7 

(1/3, -.5)  (1,0)  (-1/3,0) 

8 

(1/3, -.5)  (1.2-14)  (1,0) 

Table  1:  Example  Integration  Region  D 


For  comparison  we  include  the  results  from  TWODQ,  a single  precision 
routine  that  can  also  be  called  with  a union  of  triangles  [6].  Also  included  are 
the  results  from  TRIEX  called  over  the  eight  triangles  successively.  TWODQ 
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Tri- 

angle 

Integral 

approx. 

Est.  abs. 

error 

No.  of 
f-evals. 

Act.  abs. 

error 

1 

0.1060137D+01 

0.34D-07 

1292 

2 

0.1351354D+01 

0.52D-07 

1292 

3 

0.2846715D+01 

0.45D-06 

20872 

4 

0.1322901D+00 

0.10D-06 

4496 

TRIEX 

5 

0.3393981D+00 

0.17D-07 

224 

6 

0.1351640D+06 

0.81D-01 

402 

7 

0.1541128D+ 01 

0.48D-06 

1826 

8 

0.2427416D+03 

0.23D-03 

10548 

Total 

0.1354140D+06 

0.82D-01 

40952 

0.13D-03 

TWODQ 

ICLOSE=0 

0.1354139E+06 

0.13E+00 

40432 

0.13E+00 

ICLOSE=l 

0.1354139E+06 

0.14E+00 

64954 

0.12E+00 

DTRIST 

ICLOSE=0 

0.1354140D+06 

0.10D-00 

11872 

0.17D-01 

ICLOSE=l 

0.1354140D+06 

0.23D-01 

9568 

0.1  ID— 01 

ICLOSE=2 

0.1354141D+06 

0.91D-01 

8280 

0.48D-01 

Table  2:  Integration®  of  exp(  — * — y)((l  — x)2  + y2)  0 9 on  D 


“Requested  relative  error  for  each  integration  was  10  6. 


does  not  have  any  extrapolation  capabilities,  and  is  not  well  suited  for  the 
integration  offunctions  with  singularities.  A relative  error  of  10~* * * * * 6 *  was  requested 

for  each  integration. 

Note  that  the  triangles  1,  2,  3,  7 and  8 all  have  the  point  (1,0)  as  a vertex, 

where  the  first  integrand  has  a singularity.  In  a polar  coordinate  system  centered 

at  (1,0),  the  function  has  the  singular  behaviour  of  1 /R1'8,  where  R represents 

the  radial  coordinate  in  that  system. 

The  triangles  1 and  7 share  an  edge  along  y = 0 and  triangles  2,  3,  5,  6 and 

8 have  a vertex  on  y = 0,  where  the  second  function  has  a derivative  singularity. 
Furthermore,  the  second  function  has  higher  order  derivative  singularities  along 
* = 0.5,  which  is  inside  the  triangles  1,  2,  4,  7 and  8.  This  was  included  as  an 
example  for  which  the  extrapolation  is  not  expected  to  be  effective. 

The  TWODQ  examples  were  run  on  a Control  Data  Cyber  180/855  system 
operating  under  NOS/VE  (level  1.5.1)  at  the  National  Institute  of  Standards 
and  Technology.  This  machine  carries  about  14  decimal  digits  in  single  precision. 
The  results  from  TRIEX  and  DTRIST  were  obtained  on  a Sun  4/260  under  Unix 
4.1  at  Western  Michigan  University.  About  16  decimals  are  carried  in  double 
precision  on  this  system. 
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Tri- 

Integral 

Est.  abs. 

No.  of 

Act.  abs. 

angle 

approx. 

error 

f-evals. 

error 

1 

0.2672767D— 01 

0.12D-07 

25322 

2 

-0.4512156D  — 02 

0.36D-08 

8412 

3 

-0.1950725D^02 

0.12D-04 

224 

4 

0.5018726D+00 

0.35D-06 

6632 

TRIEX 

5 

0.1637102D+01 

0.15D-05 

402 

6 

0.1122360D+03 

0.11D-03 

402 

7 

0.131541  ID— 01 

0.33D-08 

9302 

8 

-0.2547351D+01 

0.23D-05 

3784 

Total 

0.9235573D+02 

0.12D-03 

54480 

0.19D-04 

TWODQa 

ICLOSE=0 

0.9235634E+02 

0.83E-03 

11368 

0.59E-03 

ICLOSE=l 

0.9235640E+02 

0.12E-02 

8084 

0.64E-03 

DTRIST 

ICLOSE=0 

0.9235575D+02 

0.17D-04 

6944 

0.42D-05 

ICLOSE=l 

0.9235576D+02 

0.92D-04 

9752 

0.44D-05 

ICLOSE=2 

0.9235575D+02 

0.88D-04 

12328 

0.98D-06 

Table  3:  Integration4  of  |y(0.5  — z)|1/,2(0.5  — x)  on  D 


“Requested  accuracy  was  not  reached  because  TWODQ  detected  roundoff  error. 
^Requested  relative  error  for  each  integration  was  10-6. 


5.  Concluding  remarks 

Automatic  integration  over  triangles  has  received  considerable  attention  by 
several  authors  [5, 9, 1,2, 6].  Applications  include  integration  over  irregular  re- 
gions which  can  be  triangularized  or  where  the  integral  can  be  approximated  by 
that  over  a set  of  triangles. 

TRISET  is  an  extension  of  the  TRIEX  algorithm  for  integration  over  a 
triangularized  region.  It  allows  dealing  efficiently  with  problems  where  a relative 
accuracy  is  requested  for  the  integral  over  the  entire  set  and  also  handles  some 
forms  of  boundary  singularity. 

Furthermore  TRISET  lends  itself  naturally  to  parallelization  using  the  task 
partitioning  concept  [8],  as  the  user-supplied  triangle  set  and  its  partitions  pro- 
vide an  appropriate  task  pool,  thereby  retaining  the  advantages  offered  by  the 
extrapolation  technique. 


r 
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